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We present a detailed analysis and discuss consequences of the strong correlations of the configurational 
parts of pressure and energy in their equilibrium fluctuations at fixed volume reported for simulations of sev- 
eral liquids in the previous (companion) paper |arXiv:0807.0550|. The analysis concentrates specifically on the 
single-component Lennard- Jones system. We demonstrate that the potential may be replaced, at fixed volume, 
by an effective power-law, but not simply because only short distance encounters dominate the fluctuations. 
Indeed, contributions to the fluctuations are associated with the whole first peak of the radial distribution func- 
tion, as we demonstrate by an eigenvector analysis of the spatially resolved covariance matrix. The reason the 
effective power-law works so well depends crucially on going beyond single-pair effects and on the constraint 
of fixed volume. In particular, a better approximation to the potential includes a linear term, which contributes 
to the mean values of potential energy and virial, but not to their fluctuations, for density fluctuations which con- 
serve volume. We also study in considerable detail the zero temperature limit of the (classical) crystalline phase, 
where the correlation coefficient becomes very close, but not equal, to unity, in more than one dimension; in 
one dimension the limiting value is exactly unity. In the second half of the paper we consider four consequences 
of strong pressure-energy correlations: (1) analyzing experimental data for supercritical Argon we find 96% 
correlation; (2) we discuss the particular significance acquired by the correlations for viscous van der Waals 
liquids approaching the glass transition: For strongly correlating viscous liquids knowledge of just one of the 
eight frequency-dependent thermoviscoelastic response functions basically implies knowledge of them all; (3) 
we re-interpret aging simulations of ortho-terphenyl carried out by Mossa et al. in 2002, showing their conclu- 
sions follow from the strongly correlating property; and (4) we briefly discuss the presence of the correlations 
(after appropriate time-averaging) in model biomembranes, showing that significant correlations may be present 
even in quite complex systems. 



I. INTRODUCTION 



In the companion papei^ to this work, referred to as Pa- 
per I, we detailed the existence of a strong correlation be- 
tween the configurational parts of pressure and energy in sev- 
eral model liquids. Recall that (instantaneous) pressure p and 
energy E have contributions both from particle momenta and 
positionsP 

p = NkeTip,, pn)/V + W{r,, . . . , rN)/V (1) 
£; = X(pi,...,Pw) + C/(ri,...,rjv), (2) 

where K and U are the kinetic and potential energies, re- 
spectively, and r(pi, . . . ,p7v) is the "kinetic temperature",^ 
proportional to the kinetic energy per particle. The con- 
figurational contribution to pressure is the virial W, which 
for a translationally invariant potential energy function U is 
defined-^ by 



W = - 



r,: • V.,.C/ 



(3) 



where is the position of the zth particle. Note that W has 
dimension energy. In the case of a pair potential C/pair = 
J2i<j '^{''"ij) expression for the virial becomes^ 



Impair = H njv'inj) ^-^Yl ^(""^j) ^"^^ 

i<j i<j 



where w{r) = rv'{r). 

It is the correlation between U and W that we are interested 
in, quantified by the correlation coefficient 



R 



(AWAU) 



(5) 



Paper I documented the correlation in many systems, show- 
ing that this is often quite strong, with correlation coeffi- 
cient R > 0.9, while in some other cases it was found to 
be weak or almost non-existent. The latter included in par- 
ticular models with additional significant Coulombic interac- 
tions. The purpose of this paper is two-fold. First we give 
a comprehensive analysis of the source of the correlation in 
the simplest "strongly correlating" model liquid, the standard 
single-component Lennard- Jones (SCLJ) fluid. Paper I pre- 
sented briefly an explanation in terms of an effective inverse 
power-law potential. Here we elaborate on that in greater de- 
tail, and go beyond it. Secondly we discuss a few observ- 
able consequences and applications of the strong correlations. 
These range from their measurement in a real system to ap- 
plications to systems as diverse as supercooled liquids and 
biomembranes. 

In section [11] we present a detailed analysis for the SCLJ 
case, first in terms of an effective inverse power-law with ex- 
ponent ~ 18. This accounts for the correlation at the level of 
individual pair-encounters by assuming that only the repulsive 
part of the potential, corresponding to distances less than the 
minimum of the potential r„i, is relevant for fluctuations, and 
that this may be well approximated by an inverse power-law. 
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The value 18 is significant since this explains the "slope" 7 
defined as 



(6) 



observed to be ~ 6 for Lennard- Jones systems (Paper I). The 
slope is exactly n/3 for a pure inverse power-law potential 
with exponent n, a case with perfect W, U correlation (Pa- 
per I). Section |ll] continues with a discussion of the SCLJ 
crystal, which is also strongly correlating. This would seem 
to invalidate the dominance of the repulsive part, since only 
presumably distances around and beyond the potential mini- 
mum are important, at least at low and moderate pressure. In 
this case the correlation can be explained only when summa- 
tion over all pairs is taken into account, thus the correlation 
emerges as a collective effect. There is a connection between 
the slope obtained in this way and that given by the effective 
inverse power-law; in fact they are quite similar. The third 
subsection in section [II] gives a more systematic analysis of 
which regions dominate the fluctuations in the liquid phase, 
using an eigenvector decomposition of the spatially-resolved 
covariance matrix. This matrix represents the contributions to 
the (co-)variances of potential energy and virial from differ- 
ent pair separations. It is demonstrated that the region around 
the minimum of the potential plays a substantial role. The fi- 
nal subsection in section|ll]provides a synthesis of the insights 
from the previous subsections, resulting in an "extended effec- 
tive power-law approximation", which includes a linear term. 
The main point is that a linear term in the pair potential will 
contribute to the mean values, but not to fluctuations, of W 
and U, if the volume is constant. 



In section III we discuss some consequences, starting by 



considering whether the instantaneous correlations can be re- 
lated to a measurable quantity in the normal liquid state, and 
demonstrating this with data for supercritical Argon, finding 
a correlation coefficient of 0.96. Next we focus on conse- 
quences for highly viscous liquids, where time-scale separa- 
tion implies that instantaneous correlation between virial and 
potential energy can be related to a correlation between the 
time-averaged pressure and energy. The third subsection dis- 
cusses consequences for aging, while the fourth briefly dis- 
cusses connections with recent work by Heimburg and Jack- 
son on biomembranes. 

Finally, section |lV]concludes with an outlook reflecting on 
the broader significance of strong correlation, and its implica- 
tions for the understanding of liquids, particularly in the con- 
text of viscous liquids (which has been our main motivation 
throughout this work). 



II. ANALYSIS 

A. The effective inverse power law 

In this section we consider the SCLJ system, in the hope 
that it is simple enough that a fairly complete understanding 



of the cause of strong correlations is possible. Recall that 
R > 0.9 for a wide range of states (Paper I). In order to 
understand the numerology better we consider a generalized 
Lennard- Jones potential, denoted LJ(a,b) (a > b): 



a,b 



(r) oc {a/rY - {a/rf 



(7) 



although the standard LJ(12,6) case will be used for most ex- 
amples. Starting from the idea that short distances dominate 
fluctuations and that the observed correlations are suggestive 
of a power-law interaction, we show that at a given density, 
the Lennard- Jones potential may be approximated by a sin- 
gle effective inverse power law over a range from a little less 
than a (where v^j changes sign) to around r„i, the location of 
the potential minimum (r„, = 21/6(7 for LJ(12,6)), covermg 
an energy range of approximately — e to +2e, where e is the 
well-depth. At first sight, one might expect that if this were at 
all the case, the effective power would be less than a, some- 
how a mixture of the two exponents a and b. It was noticed 
by Ben-Amotz and Stell,^ however, that the repulsive core of 
the LJ(12,6) potential may approximated by an inverse power 
law with an exponent ^18, in agreement with our data (Paper 
I). To see how we get an exponent greater than a, note that 
the exponent n in an inverse power law can be extracted from 
different ratios of derivatives: 



where _B is a constant. This implies that 



(8) 



{vpL{r)-B) 



(r) 



-'PL 



(r) 



- 2 = 



(9) 



since successive differentiation gives factors of n, then n + 
1, etc. For a potential v{r) which is not an inverse power 
law, these expressions provide different definitions of a local 
effective power-law exponent (assuming v{r) ^ as r ^ 
00): 



i(0)(r) = 



—rv'{r) /v{r) 
-rv"{r)/v'{r) - 1 
-rv"'{r)/v"{r)-2 

{r)/v^P\r) -p. 



_^y(p+i) 



(10) 



A plot of the first four of these is shown in Fig. [T] for 
LJ(12,6). All converge to a = 12 at short r, as they must. 
All increase with increasing r until the denominator vanishes, 
but more slowly, the higher the order p. In particular, when 
n'^P^ diverges it is straightforward algebra to show that n'^P+'^'i 
has the value a+b + p. So n'^^ diverges atr = a where the 
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FIG. 1: Effective power-law exponents defined by derivative ratios 
of different orders (Eqs. (|10^), for the standard Lennard- Jones poten- 
tial LJ(12,6). All converge to 12 at small r; they diverge when the 
derivative in the denominator vanishes, which happens for larger r, 
the higher the order of this derivative. The term "effective inverse 
power law" in this paper refers to a power law chosen to match "nS^^ 
at some point ro < Tm. ~ \.V2a, the potential minimum where 
•nP^^ diverges. A conve nie nt ch oice is to match at r = ct, giving 



18. In subsections 



II B 



and 



lie 



we show that n'^' (r) plays an im- 
portant role in the understanding of fluctuations associated with pair 
distances close to the potential minimum (r^ ^ 1.12(t). 



potential is zero, and is therefore unsuitable for characterizing 
the range which we expect to dominate the fluctuations, from 
energy +e to energy —t (it is also sensitive to the presence of 
an additive constant, unlike the others). Instead we use n*^^', 
which at r = (7 (the zero of v and the divergence of n'^"^ ) 
takes the value a + 6, or 18 for the LJ(12,6). Taking the factor 
3 in the denominator of the virial into account, this would ex- 
plain the slope ~6 observed in the simulations (Paper 1). But 
first we should see how well an inverse power law with this 
exponent actually fits the Lennard-Jones potential. We denote 
the point matching point rg. With the exponent fixed, we are 
free to choose the multiplicative constant A and the additive 
one B to match the slope and value at r = cr = ro; the re- 
sulting expression is (4/3)((r/cr)^^'* — 1). This is plotted in 
Fig. |2|a) along with v^j. We can match at different values 
of r by finding the expression for TiS^^(r) for the generalized 
Lennard-Jones potential: 



1 



a — h 
{b/a){r/aY-^ 



(11) 



which becomes 6 + 12/(2 - {r/af) for LJ( 12,6).™ 

The fact that we can choose a function (an inverse power 
law in this case) to match a given function and its first two 
derivatives is nothing special by itself; after all, a Taylor se- 
ries does the same. Examples of matching power laws and 
Taylor series up to fourth order, at different values of rg, 
are shown in Fig. |2|b), where the errors VLj{r) — vpL{r) 
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FIG. 2: (Color online) (a) The Lennard-Jones potential VLjir) fitted 
by an effective power law potential vpLij-) = Ar~" + B covering 
the most important part of the repulsive part of the potential. The ex- 
ponent n was chosen to be 18 which optimizes agreement at ro = cr, 
where the effective power law exactly matches not just vlj, but also 
its first two derivatives. Also shown are the Taylor series expansions 
of VLj{r) about r = up to third and fourth order. The radial 
distribution function g{r) — 1 (at T = 80K, p = 0) is also shown as 
a convenient reference for thinking about where contributions to po- 
tential energy and virial fluctuations come from. 0(b) Error made in 
approximating VLj{r) with different effective power laws matched 
at different points ro and with Taylor expansions up to third order 
about the same point. 
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FIG. 3: (Color online) Scatter-plot of true and "reconstructed" 
potential-energy and virial fluctuations (dimensionless units) for the 
LJ-liquid, where the reconstructed values UpL and WpL were cal- 
culated from the true configurations, assuming an inverse power-law 
potential with exponent 19.2; mean values have been subtracted off. 
The state point is the same as in Fig. 1 of paper I (zero average pres- 
sure, NVT ensemble). The correlation coefficients are displayed in 
the figures; the dashed lines indicate slope unity. The fact that ac- 
tual and reconstructed fluctuations correlate strongly, and with slopes 
close to unity, support the idea that the W, U correlation is derived 
from an effective inverse power-law potential dominating fluctua- 
tions. 

or VLjir) ~ VTayior{r) are plotted. The magnitude of the 
errors are similar in the range of r shown, but note that in 
the Taylor series it was necessary to match third derivatives 
at rg to achieve this, so the inverse power-law description 
is more compact. Moreover the power-law representation is 
much more useful when it comes to representing the fluctu- 
ations of total energy and virial, because an inverse power 
law (and therefore the eiTor) flattens out at a constant value at 
larger r, whereas the polynomial nature of the Taylor expan- 
sion means that away from the point of expansion, the error 
diverges rapidly (Fig.|2|a)). 

We can test the validity of the power-law approximation for 
representing fluctuations in W and U as follows. For the pur- 
pose of computing the energy and virial of a configuration due 
to a pair interaction, all necessary information is contained in 
the instantaneous radial distribution function (RDF)^ 




where p — N/V with N and V being the number of particles 
and the system volume, respectively. From this U and W may 
be computed as 

N r°° 

dr^TTr''g{r,t)vLj{r) (13) 

N r°° 

WLj{t)^~—p dr^Tir^g{r,t)wLj{,r). (14) 

Jo 



Now, UL,j{t) is re-written as an inverse power-law potential 
plus a difference term, ULj{t) = UpL{t) + Udiff{t), where 

N f°° 

UpL{t)^-p drAT:r^g{r,t)vpL{r) (15) 

Um{t) = dr4:TTr^g{r,t){vL.j{r) - vpL{r)) (16) 

^ Jo 

and similarly WLjit) = WpL^t) + Wdiff(t)- We do not in- 
clude the additive constant with the power-law approxima- 
tion; for practical reasons the potential function should be 
close to zero at a cut-off distance Tc (adding a constant to 
vpL would only add an overall constant to Upl)- We refer to 
UpL{t) and WpL{t) as "reconstructed" potential energy and 
virial, respectively, to emphasize that the configurations are 
drawn from a simulation using the Lennard- Jones potential, 
but these quantities are calculated using the inverse power law. 
In Fig. [3] we show a scatter plot of the true and reconstructed 
values of U and W, for the same state point {p = Q,T = 80K) 
as was shown in Figs. 1 and 2(a) of Paper I. Here the inverse 
power-law exponent was chosen to minimize the variances of 
the difference quantities, ((AJ/diff)^) and ((AWdiff)^)- These 
are minimized for n = 19.3 and 19.1, respectively, so we 
choose the value 19.2, which coiTesponds to matching the po- 
tentials at the distance l.OlScr. Note that what are actually 
plotted are the deviations from the respective mean values — 
the means (Ulj) and (Upl), for example, do not coincide. 
But the fluctuations are clearly highly correlated, and the data 
lie quite close to the blue dashed lines, marking slope unity. 
Specifically, the correlation coefficients are 0.946 for the po- 
tential energy and 0.984 for the virial. We can also check 
how much of the variance of Ulj is accounted for by UpL, 
((AJ7pl)^>/((AC/lj)2), and similarly for W. This is a sen- 
sible quantity because the 'PL' and 'diff parts are almost un- 
correlated for the choice n = 19.2 (cross terms account for 
less than 1% of the total variance in each case). We find 92% 
for U and 95% for W. Thus we see that the power law gives 
to a quite good approximation the fluctuations of W and U. 
The coiTelation follows from this with 7 given by one third of 
the effective inverse power-law exponent, or 6.4 for this state 
point. The measured slope (Eq. (j6])) was 6.3 corresponding to 
an effective exponent 18.9, about 2% smaller than the 19.3. A 
simpler way to determine the exponent would be three times 
the slope, although for some applications it could be advanta- 
geous to optimize the fit as described here. 

B. Low temperature limit: anharmonic vibrations of a crystal 

We turn our attention now to the fact that the correlation 
persists even for the crystallized samples (seen in Paper I in 
the lower left part of Fig. 4 and in Fig. 6). This is not trivial, 
because the physics of solids, both crystalline and amorphous 
systems, is generally dominated by fluctuations about me- 
chanically stable structures, and therefore presumably (except 
perhaps at very high pressure) by the form of the potential near 
its minimum r^, i.e., including distances larger than the min- 
imum. Thus the idea of the effective inverse power law would 
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seem to be inappropriate here — in particular since the effec- 
tive exponent n'^^ diverges at — and there is apparently no 
reason why one should get a correlation as strong as in the 
liquid and with so similar a slope. In fact there is an interval 
of r between r,„ and the minimum of the pair virial —w{r) /3 
where v{r) is increasing and —w{r)/3 is decreasing, which 
would lead if anything to a negative correlation between W 
and U, when considering individual pair interactions. More- 
over, one would expect that a harmonic approximation of the 
potential near the ground-state configuration would be an ac- 
curate representation of the dynamics in the low-temperature 
limit, but as we will see, the harmonic approximation actually 
implies negative W, U correlation, which is not observed. In 
this subsection we show why the strong positive correlation 
persists, and why the slope 7 changes little going from liquid 
to crystal (at constant volume). Although the classical dynam- 
ics of a crystal is apparently of little importance, since in real- 
ity quantum effects dominate, it turns out to be very instructive 
to consider the low-temperature {T 0) classical limit, since 
what we find has significance also in the liquid phase (subsec- 
tion |irc|. The key ideas are (1) that the positive correlation 
emerges only after summing over all interactions — it is there- 
fore a collective effect, rather than a single-pair effect, and (2) 
the constraint of fixed volume — it is important to recall from 
Paper I that the virial-potential-energy correlation only ap- 
pears under fixed-volume conditions; different volumes give 
approximately the same slope but different offsets (Fig. 4 in 
Paper I). 



1. The one-dimensional crystal 



w 



{r^.^+l - Traf] (19) 
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I 



which, by writing r^ ^+i = r,,,, + [ri^i+i — r„i), can be rewrit- 
ten as 



W=~-[k2 rraSi 



k2S2 



+ -53 



(20) 



Note that U involves 52 and 53 while W also has a first-order 
term with Si. Evaluating the sum 5i is very simple: ri.i+i can 
be expressed in terms of displacements from the equilibrium 
positions m as n^i+i = r„i + m+i - Ui, giving for 5i 



5i = ^(ui+i - Ui). 



(21) 



Such a sum of consecutive relative displacements gives the 
change in separation of the two end particles. But the total 
sum must vanish, because by periodic boundary conditions the 
"end particles" are the same particle (it doesn't matter which 
one), therefore 5i — 0. In fact periodic boundary conditions 
are not necessary, only that the length is fixed. Since both 
U and W involve at lowest order 52, which is positive semi- 
definite, at sufficiently low temperature we may drop the 53 



terms. Combining Eqs. ( 17 1 and (20i we find 



For maximum clarity we start by considering the simplest 
possible case, a one-dimensional (ID) crystal with periodic 
boundary conditions and only nearest-neighbor interactions. 
We also suppose that the lattice spacing is equal to the min- 
imum of the potential; this assumption is not made in the sub- 
sequent treatment of the 3D crystal. In a crystal the particles 
stay close to their equilibrium positions. It therefore makes 
sense to expand the pair energy (we have in mind a general 
pair potential with a single minimum) as a Taylor series, leav- 
ing out constant terms but keeping third order terms: 



^fe52 + ^^353 



(17) 



where r^.i+i is the distance between particles i and i + 1, fcp 
is the p-th derivative of the pair potential at r ~ r™, and we 
introduce the notation 



W 



l k2 + k-jTm/'i 
"3 k^ 



u 



(22) 



where we have written the coefficient in terms of the p = 2 
effective power-law exponent defined in Eq. ( [TQ| . For LJ(a,b) 
the coefficient evaluates to (a + & + l)/3, which is 6.33 for 
LJ(12,6), similar to the observed slope. This short calcula- 
tion demonstrates the main point: summing over all interac- 
tions makes the first-order term in the virial vanish, and the 
second-order term is proportional to the second-order term in 
the potential energy. It is also worth noting that for a purely 
harmonic crystal we can take fes = 0, in which case Eq. ( 22 1 
implies that there is perfect negative correlation, with a slope 
of -2/3. 



2. The three-dimensional crystal 



We now generalize this to three-dimensional crystals, 
= ^^(^i i+1 — rmY- (18) which means allowing for transverse displacements. The cal- 

j culation involves breaking overall sums into sums over one- 

dimensional chains within the crystal. We also relax the con- 
The virial is dition that the lattice constant coincides with the potential 
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minimum, which is only realistic at low pressures. We still as- 
sume only nearest-neighbor interactions are relevant (this will 
be justified in the next subsection). Generalization to a dis- 
ordered (amorphous) solid'' should be possible, since we ob- 
serve the correlation to hold also in that case. The calculation 
would necessitate, however, some kind of disorder averaging, 
which is beyond the scope of this paper.'^ 

We start by considering a simple cubic (SC) crystal of lat- 
tice constant Oc, with interactions only between nearest neigh- 
bors, so that the equilibrium bond length is Uc for all bonds. 
The fact that such a crystal is mechanically unstable is irrele- 
vant for the calculation. We shall see later that the result ap- 
plies also to, for instance, an FCC crystal. We have the same 
kind of expansions about r = Ac as above for U and W, ex- 
cept a linear term is now included, since we no longer assume 
that ac — Tm- An index h is used to represent nearest-neighbor 
bonds and as for the ID case we define 



(23) 



We then have for U and W 



p=l bonds 6 p=l 

k. 



bonds b 



-M.5o-f:(^ + ^)^, (25) 



where kp is the p-th derivative of the pair potential at r = 
flc. It is convenient to define coefficients CH and C}^ of the 
dimensionless quantities Sp/a^: 



CI 



kpaP fcp+iaj?+i 



(26) 
(27) 



Dropping the constant term —kiGcSo in W, since it plays no 
role for the fluctuations, we then have 



U^p 

p 



3H^ = ;^C7§ 



(28) 
(29) 



p=i 



The ratio of corresponding coefficients is given by the pth or- 
der effective inverse power-law exponent: 



C^/C^ =-[p+ ] = n('\a,)- (30) 



Thus for in the limit of small Cc, where these are all similar 
and close to the repulsive exponent in the potential (Fig.[T]i, the 
two expansions will be proportional to each other, to infinite 
order. Also worth noting for later use is that the CpS in each 
series increase with p. For example 



koa, 



2"c 



1 2fci 

k2ac + k^al/2 



ki + k2ac 



(31) 



(32) 



For flc between a and r„i, for LJ(12,6), the absolute values of 
these ratios lie in the intervals 9.5-cx) and 6.7-9.5 respectively. 

From dimensional considerations the variance of 5*^ is pro- 
portional to NafP, where ct^ cx T is the variance of single- 
particle displacements. At low T, therefore, we expect the Si 
terms to dominate, which causes a problem since fci changes 
sign at r„, corresponding to the divergence of n(i). In ID 
this was avoided by the exact vanishing of 5*1. In 3D 5i does 
not vanish exactly, but retains terms second order in displace- 
ments, and so contributes similarly to 82- It turns out (see 
below) that Si and 82/0.0 have similar variances and signif- 



icant positive correlation, but in view of Eqs. (31 1 and (32i 
this is not even necessary for a strong W, U correlation — the 
coefficients of the Si are relatively small so that it is still the 
^2 terms that dominate. That is essentially the explanation of 
the strong correlations in the crystal, but we now continue the 
analysis in more detail in order to investigate how good it be- 
comes in the limit T ^ 0. These general considerations will 
be of use again in the following subsection, where we make a 
similar expansion of the fluctuations in the liquid state. 

We need to evaluate 5*1 and S2 in terms of the relative dis- 
placements Ub of the pair of particles involved in bond b^^ We 
keep only terms up to second order in displacements, since we 
are interested in the limit of low-temperatures, so all ^3 terms 
in the expansion are dropped. In an SC crystal, all nearest 
neighbor bonds are parallel to one of the coordinate (crys- 
tal) axes. Consider all bonds along the x-axis. These may 
be grouped into rows of collinear bonds. The sum along a 
single row is almost analogous to the one-dimensional case, 
except that the bond length now also involves transverse 
displacements; 



fa, row 



in - a,)P 



(33) 



We write ri, explicitly in terms of relative displacements and 
expand the resulting square root, dropping terms of higher or- 
der than the second in u: 
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Similarly, for the virial 



1 H H TT- H 



1/2 




Mb,: 



(34) 



Now, the sum over bonds in a given row of the parallel dis- 
placements Mb J. vanishes for the same reasons as in the one- 
dimensional case. But when we sum the contributions to Si 
over the row, there are also second-order terms coming from 
the transverse displacements. Extending the sum to all bonds 
parallel to the x-axis, we have part of Si, denoted S'f: 



SI 



E 

b.x 



2an 



E 



2ar 



(35) 



-K2acOi — 02 



kiSi - k2S2 



kl_ 

2ar 



2 



E 1"^.- 



Ei"^ 



(40) 



3. Statistics of Si and S2 

It is clear that the || and _L sums are not equal, although 
they must be correlated to some extent; If written in terms 
of single-particle displacements rather than relative displace- 
ments for bonds, a term such as |ubj| ^ for a bond in the x- 
direction becomes 



(41) 

where R is a lattice vector used to index particles. Summing 
over bonds gives 



where _L indicates the component of the relative displacement 
vector perpendicular to the bond direction. Writing it this way 
allows us to easily include the bonds parallel to the y- and z- 
axes, and the total Si is given by 



2ar 



(36) 



S2 ~ 2 |urP — 2 (para, cross terms) 



(42) 



R 



where the cross terms are products of the parallel components 
of displacements on neighboring particles. For 5*1 we have 
something similar: 



Next we calculate 52 to second order in relative displacements 
Ub. Starting with , the part containing only bonds in the x- 
direction, using Eq. ([34| we get 



(37) 



where || means the part of the relative displacement that is 
parallel to the bond. Including all bonds. 



'^2 = ^ |Ub,| 



(38) 



Now we can return to our expressions for the potential- 
energy fluctuations (Eqs. (|24]l and p5]l), keeping only terms 
in 5*1 and 5*2: 



AU = kiSi + ^S2 = ^Y.\^tj 
2 2ar — ' 



E 



(39) 



S'lflc = 2 ^ |urP — 2 ^ (perp. cross terms) (43) 

R 



where here the cross terms involve transverse components. 
Since the term J^n I^Rp appears in both Si and ^2, these are 
correlated to some extent, but not 100% since different cross 
terms appear (note that if it were 100%, then W and U would 
both be proportional to 5*2 oc 5*1 and also correlated 100%). 
Before considering to what extent they are correlated, let us 
see how much of a difference it makes. Suppose the quanti- 
ties 5*1 Oc and ^2 have variances and ct|, respectively, and 
are correlated with correlation coefficient Rs- Using the co- 



efficients introduced in Eqs. E6s and (27 1 



Ual = CYiSiac) + C^S2 



3Wai = C^iSiQc) 



O2 >J2- 



(44) 



From this we obtain an expression for the W, U correlation 
coefficient by forming the appropriate products and taking av- 
erages: 
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R 



(45) 



a 0.998 



.2 0.996 



o 0.994 



0.992 




0.9 1 1.1 

nearest neighbor distance a^ 



FIG. 4: (Color online) Plots of predicted U correlation coef- 
ficient for T ^ for a crystal of LJ(12,6) particles for different 
degrees of correlation between the quantities Siac and 52, and of 
low-temperature simulation data. The first three curves (counting 
from the bottom) assume the variances of 5*1 ac and 5*2 are equal, and 
their correlation coefficient Rs is 0, 0.5 and 0.75, respectively. The 
fourth curve (up-triangles) results from considering a simple cubic 
lattice and assuming individual particles have uncorrelated Gaussian- 
distributed displacements, leading to specific values for the variances 
and covariance of Sia^ and 5*2. The fifth (left- triangles) shows the 
same estimate for an FCC lattice. The right-triangles are data from an 
NVT simulation of a perfect FCC crystal at T = 0.0002K. The con- 
clusion from this figure is that R does not tend to unity as T — > 0, al- 
though it becomes extremely close. The inset shows the correspond- 
ing slopes 7 (Eq. l|6]l). 



This estimation of R is plotted in Fig.|4]as a function of lat- 
tice constant for Rs —0, 0.5 and 0.75, for the case (Ti = <T2- 
Clearly the value of Rs makes little difference in the region 
of interest, Uc ^ r„i or less, where R is above 0.99. Note that 
all curves drop dramatically as the lattice constant approaches 
the inflection point (fc2 — 0) of the potential (the precise value 
at which R becomes zero depends on the statistics of S'lCc 
and 5*2). In this regime, however, higher-order terms in dis- 
placements, including ^3, ^4, etc., become more important. 



and because of Eq. ( 30 1 their inclusion tends to restore i? to a 
high value (we have not calculated their effect in detail). Also 
plotted is the estimation of R obtained by assuming that par- 
ticle displacements are uncoiTelated and Gaussian distributed 
with variance cr^ for each (Cartesian) component, correspond- 
ing to an Einstein model of the vibrational dynamics. In this 
case tedious, but straightforward algebra allows the means and 
(co-)variances of 6*1 and S2 to be calculated explicitly for 
a given lattice. The results for SC and FCC are given in Ta- 



TABLE I: Statistics of Sia^ and S2 assuming uncorrelated particle 
displacements with variance crj for each Cartesian component, for 
simple cubic (SC) and face-centered cubic (FCC) lattices. 





SC 


FCC 






UNal 


{S2) 


6Nal 


12Nal 


var(5'i ac) 


SONat 


108iV 


var(S'2) 


36Nat 


120iV 


cov(S'iac, 5*2) 




96N 


Rs 


0.73 


0.84 



ble [l] Notice that the variance of S2 is somewhat larger than 
that of S'lflc, while their means are equal. This can be traced 
to the fact the latter contains twice as many cross terms as the 
former, and a factor of one half, so the contribution from such 
terms to the mean is the same in both cases, while the con- 
tribution to the variance is smaller for SiUc- From the (co)- 
variances we find the correlation coefficient Rs = 0.73 and 
Rs = 0.84 for the SC and FCC cases respectively. These are 
also plotted in Fig. [4] A more exact calculation would take 
into account the true normal modes of the crystal, but would 
yield little of use: Data from the crystal simulations at very 
low temperature, also plotted in Fig.|4]agree within estimated 
numerical errors with both the SC and FCC estimates. The 
key point of this figure — that R is close to but less than 1 — 
apparently would change little by taking the true crystal dy- 
namics into account. In particular it is important to note that 
if i? = 1 exactly, then this would be true no matter what kind 
of weighted average of configurations is taken (what kind of 
ensemble), so a value less than unity in the Einstein approxi- 
mation is sufficient to disprove the hypothesis that i? — > 1 as 



4. The role of the coefficients Cj 



2 



Since the detailed statistics of Si and 5*2 have little effect on 
the U correlation, it must be mainly due to the numerical 
values of the coefficients C^'^ ■ We can estimate the effect of 
these by assuming Siac and S2 have equal variance and are 



uncorrelated {Rs — 0). Then according to Eq. (45 1 the W, U 
coiTelation coefficient is 



R = 



w 
2 



(46) 



which has the form of the cosine of the angle between two 
vectors = {CY ,C^) and C^' = {0^,0^)- Thus the 
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closeness of R to unity indicates that these vectors are nearly 
parallel. The tangents of the angles these vectors make with 
the Ci axis in (Ci, C2)-space are given by Eqs. (31 1 and ([32|l; 



clearly the two angles become equal in the limit of small Oc 
where 71 '^^^ and n^^^ converge. On the other hand, for Oc ~ r„ 
where ki — and n^^'> diverges, the two vectors are 



C^ = (0,fc2/2)a^ 

^-{k2,k2 + haj2)al = 

^k.ali-l,'^). 



-M^(l,l-2(n(2)+2)) 
(47) 



Clearly C'^ is parallel to the C2 axis while deviates from 
it by an angle of order 2/n^^'> ^ 1/10. The W, U correlation 
coefficient is then R = cos(l/10) - 1 - ^(I/IO)^ - 0.995, 
in agreement with the bottom curve (circles) in the main part 
of Fig.|4] In this case (a^ = r^, fci — 0, 6*1 and S2 uncorre- 
cted with equal variance), we can obtain a simple expression 
for the slope 




consistent with the result from the one-dimensional case. 

Thus when we look at the "collective" correlations in the 
crystal we naturally get a slope involving the effective power- 
law exponent n'-^h Since the latter evaluated at the poten- 
tial minimum is similar to n'^^^ at the zero of the potential, 
the slope is similar to that seen in the liquid phase. On the 
other hand, it is more typical to think about crystal dynam- 
ics starting from a harmonic approximation, adding in anhar- 
monic terms when necessary for higher accuracy. How does 
that work here? If we set k^ — as well as fci = 0, so we 
consider the purely harmonic system with the nearest neigh- 
bor distance at the minimum of the potential, then we have 
= (0,fc2/2) and C^' = -(fc2/3)(l, 1). These are not 
close to being parallel, so the correlation will be weaker (com- 
ing mainly from that of Si and S2), but more particularly, it 
will be negative, thus qualitatively different from the anhar- 
monic case. Thus the presence of the k^ affects the results at 
arbitrarily low temperature, so the harmonic approximation is 
never good enough. This is reminiscent of thermal expansion, 
which does not occur for a purely harmonic crystal. In fact, 
the Griineisen parameter for a ID crystal with nearest neigh- 
bor interactions may be shown^ to be equal to 1 + n'^^ (ac)/2. 

Finally we consider the more realistic FCC crystal. First 
note that Eqs. (36i and (38 1 are unchanged as long as Oc is 



12 nearest neighbors, four located in each of three mutually 
orthogonal planes. Taking the xy and parallel planes first, the 
neighbors are located along the diagonal directions with re- 
spect to the cubic crystal axes. As before we can do the sum 
first over bonds forming a row, then over all parallel rows. 
For a given plane there are two orthogonal sets of rows, but 



now interpreted as the nearest-neighbor distance rather than 
the cubic lattice spacing: Each position in an FCC lattice has 



the form of the sums in Eqs. ( 36 1 and ( 38 1 includes all bonds. 
The results of the calculation of (co)-variances of Si and ^2 
in the Einstein model of the dynamics are changed in a way 
that in fact increases their mutual correlation and therefore the 
W, U correlation, as shown in Table|l]and Fig.|4] 

To summarize this subsection, the correlation in the crystal 
is an anharmonic effect that persists in the limit T — > 0. It 
works because (1) the constraint of fixed volume causes the 
terms in U and W that are first order in particle displacements 
to cancel and (2) the coefficients of the "transverse" second- 
order terms are small compared to those of the "parallel" ones, 
a fact which can be traced to the resemblance of the potential 
to a power law at distances shorter than the potential mini- 
mum. "Small" here means of order 1/10, which leads to over 
99% correlation because R is essentially the cosine of this 
quantity. In one dimension there are no transverse displace- 
ments and the correlation is 100% as T — > 0; in more than 
one dimension as T ^ the correlation is very close to unity, 
but never 100%. 

To gain a more complete insight into the fluctuations, we 
next present an analysis which clarifies exactly the contribu- 
tions to fluctuations from different distances, without approx- 
imations, now again with the liquid case in mind. 



C. Fluctuation modes 

In the last two subsections we considered single-pair effects 
(associated with r < r^) and collective effects (associated 
with r ^ r„i), respectively. In this section we focus on contri- 
butions from particular pair separations without keeping track 
of which actual particles are involved. We identify the contri- 
butions to U and W coming from all pairs whose separation 
lies within a fixed small interval of separations r; fluctuations 
in the number of such pairs generate fluctuations in the the 
contributions. By considering all intervals we can systemat- 
ically analyze the variances and covariances of U and W in 
terms of pair separation, which is the purpose of this section. 

The instantaneous values of U and W are given by Eq. ( [T3| ) 
and ([14]), generalized to an arbitrary pair potential. By taking a 
time (or ensemble) average we get corresponding expressions 
for (U) and (W) in terms of g{r) = {g{r, t)), the usual ther- 
mally averaged RDF. Now we consider the variances ((AC/)-^) 
and {(AW)"^), and the covariance (AUAW). Starting 
with, for example, AU{t) = 4:npN/2 drr"^ v{r)Ag{r, t), 
where Ag{r, t) = g{r, t) — g{r), averaging and taking every- 
thing except Ag{r, t) outside the average, we have 
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((AC/)2) 
((AM/)2) = 



oo 



oo 



(47rpiV/2)M dnrf dr2riv{n)v{r2){^g{n,t)^g{r2,t)) 
Jo Jo 

/>oo />oo 

(4^p7V/2)2 / dnr? / dr2rlw{n)w{r2){Ag{n,t)Ag{r2,t)) 











(47rp7V/2)2 / dnrl / dra r2w(ri)w(r2)(A.g(ri, i)Ag(r2, t))- 



(49) 
(50) 
(51) 



Clearly the quantity which contains the essential statistical in- 
formation about the fluctuations is {Ag{ri,t)Ag{r2,t)), the 
covariance of the RDF with itself. Its magnitude is inversely 
proportional to N, so that ((AJ7)^) is proportional to N, as 
it should be. The variances of U and W and their covari- 
ance are integrals of this function with different weightings. 
To make further progress, we integrate the integrands for the 
variances over AI x M "blocks" in ri . 7'2-space. This is equiv- 
alent to considering the energy, say, as the following sum of 
AI interval-energies: 



M 



(52) 



where the interval-energy Ua (t) is defined as the integral be- 
tween boundaries and ra+i'. 



N 

Ua{t) ^ 



dr 4TTr'^g{r, t)v{r). 



(53) 



The virial can be similarly represented as a sum of con- 
tributions from the same r-intervals, W{t) = J^tLi Wa{t). 
From now on we consider the primary fluctuating quantities 
to be Ua{t) and Wa{t) and seek to understand the correla- 
tion between their respective sums in terms of correlations 
between particular [/^s and W^s. In order to achieve a reason- 
able degree of spatial resolution, we do not make the intervals 
(blocks) too big, and choose an interval width of 0.05. This 
gives AI= 42 intervals: Ui is the contribution to the energy 
coming from pairs with separation in the range 0.85 to 0.9, 
marking the lower limit of non-zero RDF, while U42 refers to 
the range 2.9 to 2.95, marking the cutoff distance of the poten- 
tial used here. We shall see explicitly that only distances up 
to around 1.4 contribute significantly to the fluctuations. We 
denote deviations from the mean as, e.g., AUa ^ Ua — {Ua) 
with the angle brackets representing the time (or ensemble) 
average. 

We are interested in the covariance of the Ua& with them- 
selves (including (AUaAUb), a b), the Wa.s with them- 
selves, and the UaS with the Wa.s- These covariances are just 



what is obtained by integrating the integrands in Eqs. (49i 



(50 1, and (51 1 over the block defined by the corresponding 
intervals for ri and r2- These values are conveniently repre- 
sented using matrices A^^, , and defined as 



(A^^),b= (AL/aAf/b 



(54) 



and 



iA^^)ah={AWaAWb) 

(A^^),fc = (AUaAWb). 



(55) 



(56) 



Note that the sum of all elements of A^'^ is the energy vari 



ance, since it reproduces the double-integral of Eq. (49 1; sim- 
ilarly for the other two matrices: 



((Af/)^) = ^(A^%, 



a, 6 



(AUAW) = ^(A^^),f, 

a, 6 



(57) 
(58) 
(59) 



At this point we make one further transformation. Define 
new matrices A*^*^*, A^^*, A^^* by 



. uu* 



A^^/((A[/)^), 



, WW* — A WW 



/{{AWr), 



(60) 



(61) 



and 



^uw* ^ A^^/V((AC/)2)((AW^)2). (62) 

This is equivalent to normalizing the UaS by the standard de- 
viation ^((AC/)2) and the WaS by y/{{AW)^) respectively. 
The elements of A^^*, A^^'* and A^^* can be thought 
of as representing, in some sense, what fraction of the total 
(co-)variance is contributed by the corresponding block. Nor- 
malized in this way, the sum over all elements of A^^* and 
^ww* jg exactly unity and that for A^^* equals the corre- 
lation coefficient R: 



. uu* 



)ab 



^(A^^*),, = l 

a,b 



(63) 
(64) 
(65) 
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TABLE II: First ten eigenvalues of the super-covariance matrix (Eq. l|66[l). their fractional contributions to the three (co-)variances (Eqs. l|63[l- 
l |65[ l) and their effective slopes (Eq. ([68j), for the SCLJ liquid with parameters as in Fig. 1 of Paper I (p=34.6 mol/1, r=80K). Contributions 
from the dominant four eigenvectors are in boldface. The last three rows list sums of the third, fourth and fifth columns over, respectively, the 
dominant four, the first ten, and all 2M eigenvectors. The sum of the fifth column over all eigenvectors should be compared (see Eq. l |65[ >) to 
the R =0.939 listed in Table I of Paper I. 



index 


eigenval. 


L'^-var. contr. 


W-var. contr. 


corr. coeff. contr. 


etiective slope 


1 


1.73 


0.01 


0.01 


0.01 


5.38 


L 


1 


U.U4 


U.Uo 


U.Uj 




J 


1.11 


0.24 


0.15 


0.19 


4.98 


4 


0.87 


0.25 


0.25 


0.25 


6.34 


5 


0.78 


0.20 


0.14 


0.17 


5.27 


6 


0.58 


0.11 


0.17 


0.13 


7.80 


7 


0.34 


0.02 


0.05 


0.03 


10.14 


8 


0.23 


0.01 


0.03 


0.01 


13.67 


9 


0.13 


0.00 


0.01 


0.00 


116.19 


10 


0.10 


0.01 


0.00 


-0.00 


-3.63 






0.797 


0.709 


0.742 




X^l,...,10 




0.884 


0.877 


0.849 




X]l,...,2M 




1.000 


1.000 


0.938 





To make a direct analysis of all possible co-variances, we 
now construct a larger 2AI x 2M matrix by placing A^*^* and 
^ww* jjjg diagonal blocks, and A'^^* and its transpose 
on the off-diagonal blocks: 



Sup 



JJW*\T ^WW* 



(66) 



This "super-covariance" matrix contains all the information 
about the covariance of the contributions of energy and virial 
with each other. This symmetric and positive semidefinite"*^ 
matrix can be separated into additive contributions by spectral 
decomposition as 



(67) 



where Va is the normalized eigenvector whose (non-negative) 
eigenvalue is Aq — this follows from the diagonalization of the 
matrix. Thus we decompose the super-covariance into a sum 
of parts. This method of accounting for the variance of many 
variables is the basis of the technique known as Principal 
Component Analysis (PCA), which is a workhorse of multi- 
variate data analysis.^ The eigenvectors represent statistically 
independent "modes of fluctuation"; the corresponding eigen- 
value is the part of the variance within the 2M-dimensional 
space accounted for by the mode. PCA is most effective when 
the eigenvectors associated with the largest few eigenvectors 
account for most of the variance in the set of fluctuating quan- 
tities. For example in the extreme case where one eigenvector 
accounts for over 99% of the variance, we could claim that 
all the different apparently random fluctuations of the differ- 
ent contributions to energy and virial were moving in a highly 



coordinated way, such that a single parameter (say, the value 
of any one of them) would be enough to give the values of all. 

In our case we are not necessarily interested in the modes 
with the largest eigenvalues: a large eigenvalue could describe 
a fluctuation mode where the individual J7aS and Wa^ change 
a lot, but their respective sums do not; this would correspond 
to the contributions from one interval increasing while those 
in others decrease, in such a way that the total is roughly con- 
stant. What we really want are those modes which contribute a 
lot to the variance of energy and virial, and to their covariance. 
This is easy to do by summing all elements in the appropriate 
block of the matrix XaVaV^ where Va, Aq, are the normalized 
eigenvector and eigenvalue in question.EHl In Table |ll] we list 
the first ten eigenvalues in decreasing order, along with their 
contributions to the normalized variances of energy, virial, and 
their covariance — the normalized covariance being equal to 
Rw.u- In addition, an "effective slope" for each mode is ob- 
tained from the a* eigenvector as 



((AT4^)2)E;:M+iK)a _ 



(68) 

where the numerator gives the sum of virial contributions for 
that mode, and the denominator the sum of energy contribu- 
tions. The factor in front, which is numerically equal to the 
overall slope 7, accounts for the standard deviation that we 
normalized the J/qS and Wa^ by to define the matrices ^vu* ^ 
A^^* and A^'^*. 

In the above equation, it looks like 7q, is determined by the 
overall 7, whereas we could expect more the opposite, that the 
overall slope is somehow an average of the individual effec- 
tive mode slopes. It looks like this because of the normaliza- 
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tion choice we made in determining the decomposition. We 
can relate the 7q to the 7 in a more meaningful way by writ- 
ing the sums in Eqs. (63 1 and (64i in terms of the spectral 



decomposition Eq. ( 67 1: 



Eg Aa(7„/7)^(Ea<MK)a)^ 
E/3-^/3(E6<Af(^/3)fc)^ 



(69) 



(70) 



(71) 



where in the last step Eq. ( [68] l was used. Multiplying both 
sides by 7^^ we get an expression for the latter as a weighted 
average of the squares of the 7q,: 



7 



(72) 



where the weight of a given mode slope 7^ is (apart from nor- 
malization) Xa = Aa(Ea<M(''^")a)^' Combining the eigen- 
value and the square of the summed "energy part" of the cor- 
responding eigenvector 

Now we can notice that the third, fourth, fifth and sixth 
eigenvectors, to be referred to respectively as EV3, EV4, EV5, 
and EV6, account for most of the three (co-)variances (to- 
talling 0.80 out of 1.00, 0.71 out of 1.00 and 0.74 out of 0.94 
for variance of U, variance W , and correlation coefficient, re- 
spectively). These four eigenvectors are represented in Fig.|5] 
We observe that, as expected, most of the fluctuations are as- 
sociated with pair-separations well within the first peak of the 
RDF, which extends to nearly r = l.Gcr (see Fig. |2|a)). In 
fact not much takes place beyond r — l.Sfj. Interestingly, of 
the four, EV5, accounting for less that 20% of the variances, 
is the only one that directly fits the idea that the fluctuations 
take place at short distances, while the other modes extend out 
to r ~ 1.3(7, beyond even the inflection point of the potential 
(around 1.24cr). 

It is instructive to repeat the fluctuation mode analysis 
for a non-strongly correlating liquid, the Dzugutov liquid at 
T = 0.65. We do not show the full results here, but they can 
be summarized as follows. There are two modes which are 
concentrated at distances less than and around the first min- 
imum of the potential. These have slopes of 5.73 and 5.01 
and contribute a total of about 0.35-0.4 to the variances and 
correlation coefficient. Since the latter is 0.585 at this tem- 
perature, these modes account for most of it. There are four 
more modes which contribute more than 5% to the variances, 
but the slopes are quite different: -9.34,7.20, 28.43 and -0.67. 
These four modes all include significant contributions at dis- 
tances corresponding to the peak in v{r); clearly this extra 
peak in the potential and the associated peak in the pair-virial 
w{r) give rise to components in the fluctuations which can- 
not be related in the manner of an effective inverse power law, 
even though fluctuations occurring around the minimum can. 
As a result the overall correlation is rather weak. 
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FIG. 5: (Color online)Representations of the eigenvectors 3,4,5, and 
6 of the super-covariance matrix. Squares represent variation in Ua 
values for a mode; diamonds represent variation in Wa values. 



D. Synthesis: why the effective power-law works even at 
longer distances 

We can apply ideas similar to those used in the crystal anal- 
ysis to understand why the correlation holds even for modes 
active at separations larger than the minimum, why the slopes 
are similar to the effective power-law slopes, and why the ef- 
fective power law works as well as it does. Recall the essential 
ingredients of the crystal analysis: the importance of summing 
over all pairs, the fixed-volume constraint and the increase of 
magnitude of coefficients of the Taylor expansion with order 
These are equally valid here, but now we use them to constrain 
the allowed deviations in g{r) from its equilibrium value, in- 
stead of displacements from a fixed equilibrium configuration. 
Define the resolved pair-density p{r) by 



p{r) = {N/2)ATTr^pg{r). 



(73) 



The requirement that this integrates to the total number of 
pairs in the system, p{r)dr = N{N~1) /2, gives a global 
constraint on fluctuations of p{r): 
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FIG. 6: Intuitive picture of allowed and disallowed fluctuations in 
p(r): (a) is not allowed because it violates the global constraint 
J Ap(r) = 0; (b) satisfies the global constraint but not locality; (c) 
could correspond, for instance, to a single bond becoming shorter, 
but this is inconsistent with fixed volume (vanishing first moment — 
such a change cannot happen in isolation); (d) is allowed — it corre- 
sponds, for example, to a single particle being displaced towards one 
neighbor and away from another. Thus one bond shortens and one 
lengthens. 



Ap(r) must be local: a peak in Ap(r) at some r must be 
compensated by an opposite peak at a nearby location rief, 
rather than one far away (thus example (b) in the figure is not 
allowed) — this corresponds to a bond having length r in the 
actual configuration and length r,ef in Fret (Fig.|6](c)). Finally 
fixed volume implies that a fluctuation cannot involve any 
substantial change in the mean nearest-neighbor bond length. 
This may be expressed mathematically as the near vanishing 
of the first moment of Ap{r): 

I rAp{r)dr ^ 0. (75) 

J first peak 

Thus, if a particle is displaced towards a neighbor on one side, 
it is displaced away from a neighbor on the opposite side, thus 
the resulting fluctuation is expected to look like Fig. |6|d), 
which is characterized by having vanishing zeroth and first 
moments. Note that we restrict the integral to the first peak. 
The principle that fluctuations of Ap(r) must be local allows 
us to write a version of Eq. ( [74| similarly restricted: 

/ Ap{r)dr = (76) 

J fiist peak 

Equations ( [75] ) and ( [76| cannot be literally true, since there 
must be contributions from fluctuations at whatever cut-off 
distance is used to define the boundary of the first peak. For 
instance, there could be a fluctuation such as Fig 6(d) cen- 
tered just to the right of this cut-off, so that only the first pos- 
itive part was included in the integrals. We can ignore these 
contributions if we assume that the potential is truncated and 
shifted to zero at the boundary, as is standard in practice (al- 
though usually at larger distances). Then fluctuations right 
at the boundary do not contribute to the potential energy. The 
fact that the only contributions to the integral are at the bound- 
ary is a restatement of the locality of fluctuations.!^ 

Now we make a Taylor series expansion of v{r) around 
the maximum tm of the first peak of g{r), using U = 



/ Ap{r)dr = 0. 
Jo 



(74) 



A typical fluctuation will have peaks around the peaks of g{r), 
but only those near the first peak will significantly affect the 
potential energy and virial (subsection 11 Ci. We can assume 
that, for a dense liquid not close to a phase transition, almost 
any configuration T may be mapped to a nearby reference con- 
figuration Fief whose RDF is identical with the thermal aver- 
age g{r). 'Nearby' implies that the particle displacements re- 
lating the F and F^ef are small compared to the inter-particle 
spacing.'*^ These displacements define the deviation Ap{r) of 
p{r) from its equilibrium value. Mapping to Fi-ef gets around 
the absence of a unique equilibrium configuration as in the 
crystal case. 

Let us consider what restriction this places on Ap{r); these 
are illustrated in Fig.|6] Because the displacements are small. 



AU = 



Ap{r) v{rM) + ki{r - tm) 



first peak \ 

+ \k2{r-rMf + 



ill) 



As for the crystal kp is the pth derivative of v{r) at the expan- 
sion point {rM here), while Mp is the pth moment of Ap{r): 



Ap{r){r — TMYdr. 



(78) 



' first peak 

A similar series exists for W, with coefficients given by 
Eq. pTl): 



3Aiy = 



(79) 
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The moments play a role exactly analogous to the sums Sp in 
the analysis of the crystal. The near vanishing of Mi corre- 
sponds to that of ^1 in the crystal case, both following from 
the fixed-volume constraint; as there, it probably holds only to 
first order in particle displacements (except in one dimension 
where it is exact), but we have not tried to make a detailed 
estimate as we did with the crystal. Recalling that the extra 
contributions to the AI2 terms will be small anyway, in view 
of Eqs. (31 -321, we simply set Mi = 0, so the two series 
become (noting that Mq = also) 



■M, 

r 



p=2 



(80) 



M 



where the coefficients C^' are those defined in Eqs. (26 1 



and (27 1, but with vm replacing a^. The points made in su6 



section 11 B regarding the relation between the two series are 
equally valid here. At orders p — 2 and higher, correspond- 
ing coefficients are related by the n'^^^ (?'a/), which are always 
all above a = 12. We expect from dimensional considera- 
tions that the variance of Mp is proportional to Nw^p, where 
wpp ^ O.Srjv/ is the width of the first peak of g{r). Thus 
moments of higher order should contribute less, and there- 
fore A/2 should dominate, implying that the proportionality 
between AU and AW is essentially n^^\rM)- This is in the 
range 15-24 for ta/ in the range 1.05(t-1.15ct, giving slopes 
between 5 and 8, similar to those observed in the fluctuation 
modes. Unlike the low-T limit of the crystal, we cannot as- 
sume the fluctuations are particularly localized around tm, so 
it is not surprising that a range of slopes show up. Notice 
that we do not see arbitrarily high mode slopes correspond- 
ing to the divergence of n^^^r) at the inflection point of the 
potential. Rather, for modes centered there, the assumption 
that we can neglect higher moments of the fluctuations no 
longer holds and there is an interpolation between n'^'^^ (r) and 
n('^)(r) which is smaller (but still greater than a = 12). 

We can now understand also why the potential and virial 
fluctuations, as "reconstructed" using the effective power-law 
potential (Fig. [3]), agree so well with the true fluctuations, even 
though the fluctuation mode analysis shows that there are sig- 
nificant contributions from distances around and beyond the 
minimum, well away from the matching point r = l.OlScr. 
Fig.jTjshows the LJ(12,6) potential, the n = 19.2 power law 



(which gave the best fit in subsection 11 A 1, and their differ- 
ence, Wdiff(7') = VLj{r) — vpL^r). The latter is obviously 
very small and flat near the matching point, but grows signif- 
icantly in an approximately linear fashion at distances larger 
than r ^ l.OScr. In view of Eqs. ( 16 1 and (73 1, a fluctuation 
of Udifi can be written as 



AU, 



diff - 



V iiff {r) A p{r)dr, 



(81) 



which has the form of an inner product of functions. Van- 
ishing fluctuations of [/diff follows if either (1) Wdiff is identi- 




FIG. 7: The true potential VLj{r), the best effective power law 
vphij') (in the sense that the fluctuations in potential energy and 
virial and reproduced most faithfully), and their difference Vi\n{j'). 
Also shown are the projected versions v^j{r) and w|ff(r) where 
the constant and linear terms (determined over the interval 0.95ct 
to 1.4(t) have been subtracted off. It is the projected functions that 
should be compared in order to make a statement about the small- 
ness of Viifi{r) relative to VLj{r), since only the projected functions 
contribute to fluctuations of total potential energy. 



cally zero, or (2) it is non-zero, but orthogonal to the space 
of allowed Ap{r). Given that Udiff is not particularly small 
except close to the matching point (see Fig. [Tji, the fact that 
f/diff fluctuations are relatively small even though they are as- 
sociated with distances away from the matching point indi- 
cates that point (2) must hold approximately. Since allowed 
Ap{r) functions are those with no constant or linear term (see 
Eqs. ( [74| and ([75|), functions orthogonal to these are those 
with only a constant and linear term: f{r) = C/o(r)+Z?/i(r) 
where fo{r) is a constant function and /i(r) is a linear func- 
tion with zero mean over the range of interest. It is clear in 
Fig. [7] that t^diff is not exactly of this form, but it can be well 
approximated by such a function. This approximation can be 
checked by standard methods of the linear algebra of function 
spaces. First we choose a range (ri, over which functions 
are to be defined. For purposes this should include the range 
of significant contributions to W and U (Fig. |5]l. We choose 
ri — 0.95a and r2 = lAa. The normalized, mutually orthog- 
onal basis vectors fo{r) and /i(r) are then given by 



Mr) = ^Nr-i - Ti 
hir) 
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(ra - riY 



(r-(ri+r2)/2). 



(82) 



The part of vnsir) which is spanned by these basis func- 
tions is -wo/o + vifi where Vi = J^^ fi{r)vdis{r)dr is the in- 
ner product of Wdiff and the coiTesponding basis vector We 
define v^^ff{r) as the part of t;diff('') projected onto the space 
of allowed functions: 
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^^diff('') 



Wdiff(»-) - vofo{r) - wi/i(r). 



(83) 



This function is also plotted in Fig. [7] where it can be seen 
that it is certainly small compared to Vdiff{r) itself. More im- 
portantly, it is also small compared to the projected part of 
VLjir), v^j{f)^ defined analogously, because it is this that 
explains why the fluctuations of [/ditf are small compared to 
those of Ul J (or equivalently Upi)- This may be quantified 
by noting that the ratio of their norms is 0.09, which indicates 
how orthogonal Udiff(?') is to the space of allowed A/9(r). If we 
projected out only the constant term from fdiff(f') and 
(the a priori more obvious way to compare the size of two 
functions) the ratio of norms would be 0.50, and it would not 
be obvious why vpL does as good a job as it does. Thus, 
again, constant volume constraint, implying Eq. ( [75| ), is im- 
portant. 

The above discussion applies equally well to the virial. We 
can now write a more accurate approximate expression for 
VLj{r), which we call the "extended effective power-law ap- 
proximation": 



(84) 



where A, B and C are constants. The associated pair virial 
{■w{r) = rv'{r)) is then 



WLj{r):^—nAr + Cr, 



(85) 



which has the same form. In both cases the term Cr con- 
tributes to the mean value but not the fluctuations, because 
/ rp{r)dr is non-zero, while / rAp{r)dr ~ for those 
Ap(r) which are allowed at fixed volume. Note also the con- 
tribution to the mean values from Cr will depend on volume 
because g{r), and hence p{r), do. Thus we can see that al- 
though there are significant contributions to fluctuations away 
from the matching point where the power law fits the true po- 
tential well, these are essentially equal for both the power-law 
and the true potential, because the difference between the two 
potentials in this region is almost orthogonal to the allowed 
fluctuations in p{r). This also explains why the fluctuation 
only holds at fixed volume (which would not be explained 
by the assumption that short-distance encounters dominate the 
fluctuations). 

The extended power-law approximation, determined em- 
pirically by the projection procedure, provides an alternative 
way to understand why the effective exponent n^^-* evaluated 
at ?■ ^ (7 agrees well with n^^^ evaluated around the mini- 
mum rjn- For the extended effective power-law approximate. 



Eq. ( 84 1 we get 



,(1) 



(r) 



-n(n+ l)Ar-("+i) 
-nAr-("+i) + C 



1 



(2). ^ ^ n(n+l)(n + 2)Ar-("+2) _ 



(86) 



Note that n(^^(r) is constant, and equal to the exponent n 
of the power-law while n(^)(r) only approaches n when r is 
small enough that the C in the denominator can be neglected 
(for the true potential n'^^^(r) increases with r and eventu- 
ally diverges, see Fig. [TJ. This emphasizes the greater use- 
fulness of n^^^ in identifying the effective power-law expo- 
nent. Recall also that our analysis earlier in this section indi- 
cates that n'^^^r), involving the second and third derivatives 
of v{r) near its minimum, is more fundamentally the cause 
of the W, U correlations, explaining something like 80% of 
the correlation in the liquid phase and over 99% in the crystal 



phase. The fact that Eq. ( 84 1 is a good approximation for the 



Lennard- Jones potential pushes the correlation to over 90% 
also in the liquid phase. 

To summarize the last two subsections, we have shown here 
that the source of the fluctuations is indeed pair-separations 
within the first peak, although only a relatively small fraction 
of the variances come from the short-r region where the ap- 
proximation of the pair potential by a power law is truly valid. 
We have also seen how the Taylor-series analysis (which in- 
volves the crucial step of taking a sum over all pairs) may 
be extended to cover the whole first peak area, with all terms 
giving roughly the same effective slope, given essentially by 
the second-order effective exponent: 7 ^ n*^^^ (r„j)/3. The 
fact that this matches the first-order effective exponent at the 
shorter distance r ^ cr is equivalent to the extended effec- 
tive power-law approximation (Eq. ([84|), which given con- 
stant volume is what justifies the replacement of the potential 
by a power law (empirically demonstrated in Fig. [3]). 



III. SOME CONSEQUENCES OF STRONG 
PRESSURE-ENERGY CORRELATIONS 

This section gives examples of consequences of strong 
pressure-energy correlations. The purpose is show that these 
are important, whenever present. Clearly, more work needs to 
be done to identify and understand all consequences of strong 
pressure energy-correlations. 



A. Measurable consequences of instantaneous 
W, U -correlations 

The observation of strong W, U correlations is of limited 
interest if it can only ever be observed in simulations. How 
can we make a comparison with experiment? In general, fluc- 
tuations of dynami cal va riables are related to thermodynamic 
response functions,IMnifor example those of U are related to 
the configurational part of the specific heat, C™"^. The lat- 
ter is obtained by subtracting off the appropriate kinetic term, 
which for a monatomic system such as Argon is SNks/^. 
The virial fluctuations, however, although related to the bulk 
modulus, are not directly accessible, because of another term 
that appears in the equation, the so-called hypervirial, which 
is not a thermodynamic quantity.^ Fortunately this difficulty 
can be handled. 

Everything in this section refers to the NVT ensemble. First 



16 



we define the various response functions and configurational 
counterparts, the isothermal bulk modulus Kt, Cv, and the 
"pressure coefficient", Py: 



Kt 
Cv 

I3v 



^ (d{p)\ 

-\dT 



Cconf 
V 

^conf 



; Kt - 
Cv~ 

P 



V 

~ iNkB 

_ NkB 

V 

NksT _ W 
V ~ V ■ 



(87) 



We also define cy = Cy /V. The following fluctuation for- 
mulas are standard (see for example Ref.|2]l 



keTV 

am') 

(AUAW) 



NkeT {W) 



V 



V 



Kj 



V 



= Cv- ^NkB = C\ 



'conf 



V(3v - Nki 



(88) 
(89) 
(90) 



Here X is the so-called "hypervirial", which gives the change 
of virial upon an instantaneous volumetric scaling of posi- 
tions. It is not a thermodynamic quantity and cannot be deter- 
mined experimentally, although it is easy to compute in simu- 
lations. For a pair potential v{r), X is a pair-sum. 



(91) 



pairs 



where x{r) = rw'{r). If we use the extended effective power- 
law approximation (including the linear term) discussed in the 



last section, then from Eq. ( 85 i we get x{r) ~ n Ar " + Cr. 



Summing over all pairs, and recalling that when the volume is 
fixed the Cr term gives a constant, we have a relation between 
the total virial and total hypervirial. 




6.2 



W/N 



FIG. 8: Scatter-plot of instantaneous virial and hypervirial (in di- 
mensionless units) for a SCLJ system at p = 1.0, T = 0.80 (NVE). 
The correlation coefficient between these quantities is 0.998. The 
hypervirial is the main contribution to the configurational part of the 
bulk modulus; it gives (after dividing by volume) the change of virial 
for a given relative change in volume. The sizeable constant term in 
the linear fit shows that Eq. ^92\ is a poor approximation. The slope is 
4.9, about 10% smaller than 7 ~ 5.4 for this state point. The differ- 
ence reflects the limit of the validity of the power-law description-in 
fact a more detailed analysis shows that the relation between W and 
X is dominated by n'^^ (r), which is smaller than n^^' (r) (Fig.[Tjl. 



Inserting the fluctuation formulas Eqs. ( 88 1 and ( 90 1 gives 



R^{{p)-Kt 



V ' 



1 (keT'^VPT^f 



keTV 



keT'^C^"^ 



= T 



conf \ 2 



^conf 



(95) 
(96) 



X = {n/3)W + constant. 



(92) 



This constant survives, of course, when we take the ther- 
mal average {X), as do the corresponding constants in (U), 
(W). To get rid of these constants one possibility would be to 
take derivatives with respect to T, but this can be problematic 
when analyzing experimental data. Instead we simply com- 
pare quantities at any temperature to those at some reference 
temperature Ty^f, this effectively subtracts off the unknown 
constants. Taking first the square of the correlation coeffi- 
cient, we have 



Defining quantities A = (p) - Kt + {X) /V and B = 
r(/?™"f)2/cy°^ (the reason for the tilde on A will become 
clear), we have E?A = B. This is an exact relation. To deal 
with the hypervirial we first take differences with the quanti- 
ties at Tref, assuming that the variation of R is much smaller 
than the A and B variations: 



R^{A-A,f) = B- 



B 



'ref 



(97) 



where A,-ef = A{Ty^f), etc. A — A^^f written out explicitly is 



which implies 



((At/AVF))^ 
((Al/)2)((AVF)2)' 



R 



{{AUAW)f 



ksTV keTV ((AC/)2) 



(93) 



(94) 



A~A„ 



{{p)-KT)-{{pU-KT,,i) + 



{X) (X) 



ref 



V 



(98) 

Next we use the power-law approximation to replace {X) — 
(Xref) with {72/3){{W) - (W,ef)). This is the crucial point: 
Whereas it is often not a good approximation that {X) = 
{n/3){W) due to the unknown additive constants discussed 
above, subtracting two state points considers changes of {X) 
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and (W) with temperature. Recall from section III B of Pa- 
per 1 that the changes in mean-values A{W) and A{X) be- 
tween (nearby) temperatures are related as the linear regres- 
sion of the fluctuations of those quantities at one (nearby or 
intermediate) temperature. The linear relation between sub- 
tracted mean values holds if the instantaneous W and X are 
strongly correlated in the region of interest. The latter is con- 
firmed by our simulations; indeed the correlation of instanta- 
neous values of X and W is even stronger than for W and 
U, with approximately the same slope (Fig. [8]). Thus Eq. ( [98] l 
becomes 



A - Aef -Hp)- Kt) - ((p)ref " i^Tref) + 



n {W) - {W) 



ref 



V 

(99) 
(100) 



where A = p—KT+{n/?>){{W) /V) (no tilde) contains quan- 
tities which are all directly accessible to experiment, except 
for the effective power-law exponent n. This can be obtained 
by noting that if there were perfect correlation, one could in- 
terchange AW and (n/3)A[/; thus 



,conf 



{AUAW) 
((AC/)2) 



(101) 



which gives for A 



A^{p)-Kt + 



^p^conf^conf 
^conf 



(102) 



Thus to compare with experiment one should plot B — i?,ef 
against A—Aj^f, the prediction, in the case of near perfect cor- 
relation, ~ 1, is a straight line with slope close to unity. 
Fig. |9] shows data for Argon for T between 200K and 660K. 
Argon was chosen because as a monatomic system there are 
no rotational or vibrational modes contributing to the heat ca- 
pacity and it is therefore straightforward to subtract off the 
kinetic part. Also we restrict to relatively high temperature 
to avoid quantum effects. The correlation coefficient R given 
as the square root of the slope of a linear fit is 0.96. Note that 
the assumed constancy of R is confirmed (going to lower tem- 
peratures there are increasingly large deviations). The impor- 
tance of subtracting from a reference state point is highlighted 
by the inset, which shows that A{T) = B{T) does not hold: 
There is a correlation in the fluctuations which is not present 
in the full equation of state .l^ 



B. Time averaging: Pressure-energy correlations in liighly 
viscous liquids 

We have observed and discussed in paper 1 that when vol- 
ume is held constant, the correlations tend to become more 
perfect with increasing T, while along an isobar (considering 
still fixed-volume simulations, choosing the volume to give a 



Supercritical Argon (200K < T < 660K) 



A = <p> - + <p >(3y /Cy 



b=t(p7'')V;"' 



o 



« 0.1 



o 



0.01 



0.01 




R = 1.00 

-- R = 0.96 
O 20 mol/1 
25 mol/1 
30 mol/1 



A(700K)-A(T) [GPa] 



FIG. 9: (Color online) Data from the NIST databastP for su- 
percritical Argon at three different densities covering the temper- 
ature range 200K-660K showing a strong virial-potential energy 
correlation (i? = 0.96) [reproduced from Ref. i4j|. Here Kt = 
-V{d{p)/dV)T, p'^™' = P ~ NksT/V = W/V /^r' = 
{d{p)/dT)v - Nkg/V, and c^' = Cv/V ~ (3/2)iVfcfl/V The 
diagonal line corresponds to perfect correlation. The inset shows 
"unsubtracted" values for A and B; the fact that the data do not fall 
on the solid line indicates that a power-law description does not hold 
for the full thermodynamics. 



prescribed average pressure), they become more perfect with 
decreasing T. This fact makes the presence of correlations 
highly relevant for the physics of highly viscous liquids ap- 
proaching the glass transition. Basic questions such as the ori- 
gins of non-exponential relaxations and non-Arrhenius tem- 
perature de pendence are still vigorously debated in this field 
of researchll2Ein3IiIl Instantaneous correlations of the kind 
discussed in this work would seem to be relevant only to the 
high frequency properties of a highly viscous liquid; their rel- 
evance to the long time scales on which structural relaxation 
occurs follows from the separation of time scales as explained 
below. 

A question that is not actively debated in this research field 
(but see, e.g., Refs. 16 1718), is whether a single parameter is 
enough to describe a highly viscous liquid. The consensus for 
more than 30 years is that with few exceptions these liquids 
require more than just one parameter, a conclusion scarcely 
surprising given their complexity. The meaning of "having a 
single parameter" can be understood as follows. Following a 
sudden change in volume, both pressure and energy relax to 
their equilibrium values over a time scale of minutes or even 
hours, sufficiently close to the glass transition. If a single pa- 
rameter governs the internal relaxation of the liquid, then both 
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pressure and energy relax with the same time sca le, an d in fact 
the normalized relaxation functions are identicalP^E^This be- 
havior can be expressed in the frequency domain, as a certain 
quantity, the dynamic Prigogine-Defay ratio, being equal to 
unity. '- A key feature of highly viscous liquids is the separa- 
tion of time scales between the slow structural ("alpha") relax- 
ation (up to of order seconds) and the very short times (of or- 
der picoseconds) characterizing the vibrational motion of the 
molecules. This separation allows a more direct experimental 
consequence of W, U correlations than that described in the 
previous subsection: Suppose a highly viscous liquid has per- 
fectly correlated W, U fluctuations. When W and U are time- 
averaged over, say, one tenth of the alpha relaxation time t^T^ 
they still correlate 100%. Since the kinetic contribution to 
pressure is fast, its time-average over rQ,/10 is just its thermal 
average, and thus the time-averaged pressure equals the time- 
average of W/V plus a constant. Similarly, the time-averaged 
energy equals the time-averaged potential energy plus a con- 
stant. Thus the fluctuations of the time-averaged W and U 
equal the slowly fluctuating parts of pressure and energy, so 
these slow parts will also correlate 100% in their fluctuations. 
In this way we get from the non-observable quantities W and 
U to the observable ones E and p (we similarly averaged to 
observe the correlation in the SQW system in Paper I). The 
upper part of Fig.[TO]shows normalized fluctuations of energy 
and pressure for the commonly studied Kob-Andersen binary 
Lennard-Jones system^*' (referred to as KABLJ in Paper 1), 
time-averaged over one tenth of Tq,. In the lower part we show 
the dynamic Prigogine-Defay ratio, which in the NVT en- 
semble is defined as follows: 



and 



(103) 



here kt — is the isothermal compressibility and " 

denotes the imaginary part of the complex frequency de- 
pendent response function. A way to interpret this quantity 
can be found by considering the fluctuation-dissipation theo- 
rem expressions for the response functions. For example the 
frequency-dependent constant-volume specific heat cy(a;) is 
giverl^by 



Cy(w) = 



{AE{0)AE{t)) exp{-iujt)dt 
(104) 



where E is the total energy. Taking the imaginary part we 
have 
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where we use C to represent Laplace transformation. Simi- 
lai'ly 
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Forming the Prigogine-Defay ratio then gives, after cancelling 
factors of fc^, T, V, and uj, 



{C{{AE{0)AEmy{C{{Ap{0)Apmy 
iC{{AE{0)Apmy^ 



ATy(w) = 

\i^r,[uiLj>.p[ii/ f) - 

(108) 

We can see that the right-hand side has a similar structure to 
a correlation coefficient, if we take the inverse square root. 
So in a loose sense the dynamic Prigogine-Defay ratio can be 
thought of as the inverse square of a correlation coefficient, re- 
ferred to a particular time scale. This gives an intuitive reason 
for why it is in general greater or equal to unity, with equality 
only achieved in the case of perfect correlation.'** The lower 
panel of Fig. [T0| shows this quantity for a range of frequencies 
for the KABLJ system. It clearly approaches one at low fre- 
quencies, and stays within 20% of one in the main relaxation 
region. In the sense above, this corresponds to i? > 0.9, or 
strongly correlation. 

The line of reasoning presented here opens for a new way of 
utilizing computer simulations to understand ultraviscous liq- 
uids. Present-day computers are barely able to simulate one 
microsecond of real-time dynamics, making it difficult to pre- 
dict the behavior of liquids approaching the laboratory glass 
transition. We find that pressure-energy correlations are al- 
most independent of viscosity, however, which makes it pos- 
sible to make predictions regarding the relaxation properties 
even in the second or hour range of characteristic times. Thus 
if a glass-forming liquid at high temperatures (low viscosity) 
has very strong pressure-energy correlations (R ^ 1), its eight 
thermoviscoelastic response functions at ultraviscous condi- 
tions may basically be expressed in terms of just one, irre- 
spective of temperature (or viscosity). 



C. Aging and energy landscapes 

We now discuss the significance of the present results for 
the interesting results reported in 2002 by Mossa et fl/,'^who 
studied the inherent states visited by the Lewis-Wahnstrom 
model^^^ of the glass-forming liquid orthoterphenyl (OTP) dur- 
ing aging, i.e., the approach to equilibrium. An inherent 
state (IS) is a local minimum of the so-called potential energy 
landscape (PEL) to which a giv en co nfiguration is mapped 
by steepest-descent minimiz ation .I21|25| -j^g PEL formalism in- 
volves modeling the distributions and averages of properties 
of the IS in the hope of achieving a compact description of 
the thermodynamics of glass-forming liquids. The thesis 
of Ref. 22 and of the previous work^^ is that an equation of 
state can be derived using this formalism which is valid even 
for non-equilibrium situations. This involves including an ex- 
tra parameter, namely the average IS (potential) energy, (e/5), 
so that the equation of state takes the form 
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FIG. 10: (Color online) Upper panel, time-averaged (over Ta/10, 
where r„ is the structural relaxation time), normalized fluctuations 
of E and p in NVT simulations of the Kob- AnderseiP 

binary 

Lennard-Jones ( KABLJ) system, plotted against time in units of 
Ta ~ 10'^ (J AA\/m/eAA- The density was 1.2cr^^, and the temper- 
ature 0.474eAA- Middle panel, imaginary parts of the three response 
functions — c„(ci;), —I3v{i^) and l/KT('i'), scaled to the maximum 
value. Lower panel, dynamic Prigogine-Defay ratio for the same 
simulation. The approach towards unity at frequencies smaller than 
the loss-peak frequency (~ 1/ra) is exactly equivalent to the cor- 
relation between time-averaged quantities shown in the upper panel 
[reproduced from Ref. il9J . 



p{T,V,{eis)) = pisi{eis),V) + p,UT,V,{eis)) (109) 

where pjg is the ensemble averaged inherent state pressure — 
for a given configuration it is the pressure of the corresponding 
inherent state — andpvib = P^Pis- The usefulness of splitting 
in this way lies in the fact that pis does not explicitly depend 
on T. 

In equilibrium (e/s) is determined by V and T. The 
conclusion of Ref. |22| is that knowledge of (e/s) in non- 
equilibrium situations is enough to predict the corresponding 
pressure (given also T and V). This was based on extensive 
simulations of various aging schedules. Thus the authors con- 
clude that the inherent states visited by the system while out of 
equilibrium, must be in some sense the "same" ones sampled 
during equilibrium conditions. "Same" is effectively defined 
by their results that averages of various inherent state proper- 
ties (V, e/s, Pis, as well as a measure of the IS curvature) are 
all related to each other the same way under non-equilibrium 
conditions as under equilibrium conditions. It was similarly 
found that the volume could be determined from (e/s) fol- 
lowing a pressure jump in a pressure-controlled simulation. 
On the other hand, subsequent work by the same group found 
that this was not at all possible for glassy water during com- 
pression/decompression cycles.l22l 

Now, our results (Paper I) for the same OTP model show 
that it is a strongly coiTelating liquid. Thus we expect a gen- 
eral correlation between individual, not just average, values 
of pjs, the inherent state pressure (which lacks a kinetic term 
and therefore equals the inherent state virial divided by vol- 
ume) and e/s, for a given volume. Therefore, for any given 
collection of inherent states with the same volume — not just 
equilibrium ensembles — the mean values of U and W will fall 
on the same straight line as the instantaneous values. Note that 
this would not hold if the correlation was non-linear CoiTe- 
spondingly for given p/s, there is a general coiTelation be- 
tween individual values of e/s and V. In fact, any two of 
these quantities determine the third with high accuracy, and 
this is true at the level of individual configurations, including 
inherent states. 

To see how this works for cases involving fixed volume, we 
write the total (instantaneous) pressure as a sum of an inherent 
state part, which involves the the virial at the corresponding 
inherent state, plus a term involving the difference of the true 
virial from the inherent state virial, plus the kinetic term: 
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(110) 



The first term is linearly related to the inherent state energy for 
a strongly correlating liquid. Moreover, the difference term is 
similarly expressed in terms of the corresponding energy dif- 
ference, W — Wis = i{U — ejs)- Taking averages over the 
(possibly non-equilibiium, although we assume equilibrium 
within a given potential energy basin) ensemble, we expect 
that {U — eis) depends only on T and V (a slight e/s de- 
pendence can appear in 7 since this is slightly state-point de- 
pendent). Thus it follows that p can be reconstructed from a 
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TABLE III: Data from NpT simulations of fully hydrated phospho- 
lipid membranes of l,2-Dimyristoyl-sn-Glycero-3-Phosphocholine 
(DMPC), l,2-Dimyristoyl-sn-Glycero-3-Phospho-L-Serine with 
sodium as counter ion (DMPS-Na), hydrated DMPS (DMPSH), 
and l,2-Dipalmitoyl-sn-Glycero-3-Phosphocholine (DPPC)P^ 
The columns list temperature, correlation coefficient between 
volume and energy, average lateral area per lipid, simulation time in 
equilibrium, and total simulation time. 

T [K] Rev Aup [A^] t [ns] Uot [ns] 



DMPC 


310 


0.885 


53.1 


60 


114 


DMPC 


330 


0.806 


59.0 


50 


87 


DMPS-Na 


340 


0.835 


45.0 


22 


80 


DPPC 


325 


0.866 


67.3 


13 


180 



knowledge of (average) e/5, V and T, without any assump- 
tions about the nature of the inherent states visited. In par- 
ticular, no conclusion can be drawn regarding the latter. The 
failure of the pressure reconstruction in the case of wateiEHis 
not surprising, since water models are generally not strongly 
correlating (which as we saw in Paper I is linked to the exis- 
tence of the density maximum). 



D. Biomembranes 
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FIG. 11: (color online) Normalized fluctuations of energy (x) 
and volume (o) for a l,2-Dipalmitoyl-sn-Glycero-3-Phosphocholine 
(DPPC) membrane at 325 k!^ Each data point represents a 0.5 ns 
average. Energy and volume are correlated with a correlation coeffi- 
cient of 7? = 0.87 (NpT). 



A completely different area of relevance for the type of cor- 
relations reported here relates to recent work of Heimburg 
and Jackson,^" who have proposed a controversial new the- 
ory of nerve signal propagation. Based on experiment and 
theory they suggest that a nerve signal is not primarily electri- 
cal, but a soliton sound-wave.^' Among other things this the- 
ory explains how anaesthesia works (and why one can dope 
people with the inert gas Xenon): Anaesthesia works sim- 
ply by a freezing-point depression that changes the membrane 
phase transition temperature and affects its ability to carry 
the soliton sound wave. A crucial ingredient of the theory 
is the postulate of proportionality between volume and en- 
thalpy of microstates, i.e., that their thermal equilibrium fluc- 
tuations should correlate perfectly. This should apply even 
through the first-order membrane melting transition. The the- 
ory was justified in part from previous experiments by He- 
imburg showing proportionality between compressibility and 
specific heat through the phase transition.^- The postulated 
correlation — including the claim that is survives a first-order 
phase transition — fits precisely the pattem found in our liquid 
simulations. 

By re-examining existing simulation as well as 

carrying out extensive new simulations,'^! we have investi- 
gated whether the correlations are also found in several model 



membrane systems, five of which are listed in Table III The 



simulations involved a layer of phospholipid membrane sur- 
rounded by water, in the La phase (that is, at temperatures 
above the transition to the gel-state), at constant p and T. 
When p, rather than V, is constant, the relevant quantities that 
may correlate are energy and volume. As with viscous liq- 



uids (subsection III B 1 and the square well system (Paper I), 
time-averaging is necessary for a correlation to emerge, where 
now 



AE{t) ~ 7™'AT/(t) 



(111) 



When an averaging time of 1 ns is chosen, a significant cor- 
relation emerges, with correlations between 0.8 and 0.9 (Ta- 



ble III note these Rev values fall slightly short of our cri- 
terion of 0.9 for "strong correlations"). The time series of 
time-averaged, normalized E and V for one case are shown in 
Fig.[TT] The necessity of time-averaging stems from the pres- 
ence of water, which we know does not exhibit strong corre- 
lations. Since the membrane dynamics are much slower than 
those of the water, they can be isolated by time averaging. 



IV. CONCLUSIONS AND OUTLOOK 

In Paper I and this work we have demonstrated several cases 
of strongly correlating liquids and some cases where the cor- 
relation is weak or absent (at least under normal conditions 
of pressure and temperature). An important next step is to 
continue to document the existence or non-existence of corre- 
lations, particularly in different kinds of model systems, such 
as non-pair potential systems and systems with interactions 
computed using quantum mechanics (e.g., by density func- 
tional theory). It is noteworthy in this respect that after suit- 
able time averaging the correlations may appear in systems 
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where they are otherwise unexpected. One example was the 
square-well (SQW) case, where the correlation was between 
the time-averaged virial and potential energy. In the case of 
viscous liquids time averaging allowed a correlation to ap- 
pear between the more accessible energy and pressure, while 
for the biomembranes it made it possible to remove the non 
strongly-correlating contributions from the water In all cases 
time averaging is relevant because of a separation of time 
scales: the square-well system because the time scale over 
which the average number of neighbors changes is long com- 
pared to the time between collisions; viscous liquids because 
the vibrational dynamics (which includes the kinetic contri- 
butions) is fast compared to the slow configurational dynam- 
ics; biomembranes because the membrane dynamics is slow 
compared to those of the water Note that in this case it was 
necessary to consider an NpT ensemble and study correlations 
between energy and volume, because only a part of the system 
is strongly correlating, and this part cannot be constrained to 
a particular volume. It is worth noting that even with fixed 
volume, the correlation coefficient depends on whether the 
ensemble is NVT or NVE, although the strongly correlating 
limit of i? 1 is independent of ensemble. 

A point which has been mentioned, but which is worth em- 
phasizing again, is that the replacement of the potential by an 
appropriate inverse power law can only reproduce the fluctua- 
tions, and not the mean-values of (potential) energy and virial, 
nor their first derivatives with respect to T and V . These de- 
termine the equation of state, in particular features such as 
the van der Waals loop which are absent in a pure power- 
law system, even if changes in the exponent are allowed. The 
generalization to the extended effective inverse power-law ap- 
proximation, however, allows in principle such features to be 
described. We are currently investigating the dependence of 
the parameters of the extended approximation on state point. 

Finally, consider again viscous liquids, which are typi- 
cally deeply supercooled. The most common way of clas- 
sifying them involves the fragility parameter introduced by 
Angell,^^ related to the departure from Arrhenius behavior 
of the temperature dependence of the viscosity. Strong liq- 
uids, having the most Arrhenius behavior, have traditionally 



been considered the easiest ones to understand, because Ar- 
rhenius temperature dependence is well-understood. But it 
may well be that strongly correlating liquids, are in fact the 
simplest.^ The connection with the long-discussed question 
of whether a single-order parameter describes highly viscous 
liquids has been discussed briefly in subsection |11IB[ and is 
discussed further in Ref 16 As an example, a direct appli- 
cation of the "strongly-correlating" property concerns diffu- 
sion in supercooled liquids. Recent work by Coslovich and 
Roland^*^ has shown that the diffusion constant in viscous bi- 
nary Lennard-Jones mixtures may be fitted by an expression 
D = F{p'^ /T) where 7 reflects the effective inverse power 
of the repulsive core. "Densit y scaling" has also been ob- 
served experimentally.E2EIiIEl It is natural, given Coslovich 
and Roland's results, to hypothesize that the scaling exponent 
is connected to pressure-energy correlations, and in Ref. |4] 
it was conjectured that density scahng applies if and only if 
the liquid is strongly correlating. We have very recently stud- 
ied the relationship between the two quantitatively'*'* and have 
found that (1) density scaling does indeed hold to the extent 
that the liquids are strongly correlating, and (2) the scaling ex- 
ponent is given accurately by the slope 7 of the correlations 
(hence our use of the same symbol). This finding supports 
the conjecture that strongly correlating liquids may be sim- 
pler than liquids in general. 

In summary, the property of strong correlation between the 
equilibrium fluctuations of virial and potential energy allows 
a new way to classify liquids. It is too soon to tell how fruitful 
this will turn out in the long term, but judging from the appli- 
cations briefly presented here, it seems at least plausible that 
it will be quite useful. 
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The analysis of fluctuations in the liquid in subsection 



lie 



be applied to the case of a disordered solid, explaining the high 
correlation also there, but it is not accurate enough to get as good 
an estimate for the low-temperature limit as we do our analysis of 
the crystal. 

"Displacement" refers to the displacement of a given particle from 
its equilibrium position, while the "relative displacement" is the 
difference in this quantity for the given pair of particles. 
This follows since if v is an arbitrary vector and Aab = 
(AxaAa;;,) is the covariance matrix of a set {xa} of random 
variables, then v ■ A ■ v is the variance of the random variable 
w — "^^VaXa, and thus non-negative. 



The sums over the diagonal blocks are non-negative, since Aq is, 
and the sum over the first (last) components of Vav"^ is the dot 
product of the first (last) half of with itself, which is also non- 
negative. 

The configuration Frcf is analogous to the inherent st ate co nfigu- 
ration often used to describe viscous liquid dynamicsj^^'^ which 
is obtained by minimizing the potential energy starting from con- 
figuration F. 

We have checked the statements that the zeroth and first moments 
of p(r) over the first peak are constant apart from contributions 
at the cutoff by computing orthogonalized versions of them (us- 
ing Legendre polynomials defined on the interval (0.8(j,1.4cr)) and 
showing that they are strongly coiTelated (correlation coefficient 
0.9) with a slope corresponding to the cutoff itself. 



